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THE  GLOBAL  POSITIONING  SYSTEM 
VERSUS 

GRAVITY  DISTURBANCE  MODELING 
IN  AN 

INERTIAL  NAVIGATION  SYSTEM 


Abstract 

The  Global  Positioning  System  provides  a  real-time  means  of  updating  an  aircraft  inertial 
navigation  system  to  reduce  position  and  velocity  errors  due,  in  part,  to  an  inexact 
knowledge  of  the  gravity  field  in  which  the  aircraft  is  flying.  Gravity  disturbance 
modeling  via  such  techniques  as  point  mass  or  finite  element  modeling  provides  a 
real-time  means  of  reducing  this  gravity  field  error.  This  paper  demonstrates  that  should 
the  navigator  be  denied  the  GPS  signal,  the  modeling  of  gravity  disturbances  can 
adequately  minimize  the  navigation  error.  In  the  presence  of  the  GPS  signal,  the  gravity 
disturbance  component  information,  combined  with  GPS  data  via  Kalman  filtering, 
consl^titutes  a  further  refinement  in  the  system. 


THE  GLOBAL  POSITIONING  SYSTEM  VERSUS  GRAVITY 
DISTURBANCE  MODELING  IN  AN  INERTIAL  NAVIGATION  SYSTEM 

I  INTRODUCTION 

Recent  developments  in  inertial  instrumentation  have  substantially  improved  the 
performance  of  inertial  navigation  systems  in  both  aircraft  and  shipboard  applications. 
The  level  of  confidence  in  these  inertial  systems  has  risen  to  the  point  that  the  geodetic 
community  has  begun  to  utilize  inertiaJ  instrumentation  to  perform  surveys  of  geodetic 
quality  (l).  In  the  geodetic  survey  mode,  data  from  inertial  systems  are  processed 
post-mission  to  determine  the  components  of  the  gravity  disturbance  vector.  However, 
in  the  navigation  mode,  operating  in  real-time,  these  gravity  disturbance  components 
must  somehow  be  furnished  to  the  inertial  navigation  system,  in  which  they  have  now 
become  a  significant  error  source. 

In  state-of-the-art  inertial  navigation  systems,  electrostatically  supported  gyros 
have  reduced  drift  rates  to  as  low  as  0.00000 l°/hr  and  electrostatic  accelerometers  have 
reduced  accelerometer  bias  to  as  low  as  1  pg  (~1  x  10“  m/sec/sec)  (2).  The 
accelerometer  is  a  specific  force  instrument,  i.e.,  it  measures  the  total  acceleration  due 
to  the  vehicle,  errors  and  biases,  gravity.  Hence,  the  unknown  gravity  disturbance 
components  act  as  time  dependent  accelerometer  biases,  which,  by  Einstein's  Principle  of 
Equivalence,  cannot  be  separated  by  instrumentation  operating  at  the  acceleration  level. 

Although  the  gravity  gradiometer  shows  great  promise  for  real-time  determination 
of  the  gravity  disturbance  components  (3),  the  implementation  of  such  a  system  remains 
just  beyond  our  grasp.  Thus,  we  are  led  to  the  necessity  to  provide  either  a  means  of 
modeling  the  gravity  disturbance  in  real-time  or  a  means  of  updating  the  inertial  naviga¬ 
tion  system  with  position  and  velocity  information  in  real-time,  so  that  the  errors  induced 
by  the  unmodeled  gravity  disturbances  can  be  minimized. 


Position  information  derived  from  electronic  systems  such  as  VORTAC,  DME, 

LORAN  and  velocity  information  derived  from  Doppler  radar  is  sufficiently  accurate  for 

many  applications.  -To  satisfy  more  stringent  requirements,  the  Global  Positioning  System 

(GPS)  has  been  developed.  When  used  with  appropriate  receivers,  it  is  capable^  deriving  *  OS 

I —  /C-dP£>  j —  '—7 

a  moving  vehicle's  position  and  velocity  to  an  accuracy  o^bout  -  8  meters  /and  4  .05  I 

meters/second  (one  sigma)  (4). 

The  availability  of  the  Qobai  Positioning  System  begs  the  question  of  whether  it 
negates  the  necessity  of  modeling  the  gravity  disturbance  vector,  i.e.  --Whether  there  is 
any  need  to  model  the  gravity  disturbance,  since  we  can  cancel  its  effect  by  using  the 
GPS.  This  paper  answers  that  question  by  using  error  propagation  techniques  to  show  the 
effect  of  the  GPS  data  in  a  typical  navigation  system,  both  with  and  without  compensa¬ 
tion  by  gravity  disturbance  modeling. 

In  the  modern  inertial  navigation  system  (Figure  1)  the  gravity  vector  g  computed 
in  a  feedback  loop  contains  terms  due  to  the  earth's  gravitational  field  and  due  to  the 
earth's  rotation: 


g  =  G  -  n.  Jl.  r 
“  —  le  le- 


where 


g  =  gravity  vector 
G  =  gravitation  vector 

=  skew  symmetric  matrix  form  of  earth's  inertially  =000 
referenced  angular  velocity, 0  0 

£=  geocentric  position  vector  of  vehicle  0  oj.^  0  J 

Now,  the  gravitation  vector  G  is  itself  the  sum  of  a  reference  gravitation  vector  y 
and  the  gravity  disturbance  vector  6g; 

G  tig  (2) 

Typically,  the  inertial  navigation  system  computes  the  components  of  x  terms 
of  those  terms  which  correspond  to  the  mathematical  reference  figure  used  for  the  earth. 
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an  ellipsoid.  In  terms  of  spherical  harmonics,  the  componets  of  are  given  by: 


L-30COS  L  +  3)J 

Ye  =  0 

y  =-3GM  sinLcosL  32  ^  5  3,  R  . 

"  7~  VtA  L  ‘  ' 

(7  cos^  L  -  3)  I 


where 

GM  =  gravitational  constant  times  mass  of  Earth 
L  =  geocentric  latitude  of  vehicle  (approximated  by  geodetic  latitude) 

R  =  equatorial  radius  of  earth 

^2,34  =  second  and  fourth  degree  zonal  harmonic  coefficients  of  geopotential 

The  earth  rotation  term,  £,  is  computed  along  with  the  reference  gravita¬ 

tion  vector  so  that  we  may  consider  the  gravity  vector  as 

g=2*+l£ 

* 

where  ^  =  reference  gravity  vector. 

The  components  of  the  reference  gravity  vector  are  given  by 

y\.0  (5) 

*  2 

Y  =  Y  -  lu  •  sin  L  cos  L 
'  n  'n  le 

u).  =  earth  rotation  rate 
le 

The  gravity  disturbance  vector,  fig,  which  can  be  defined  via  equation  (4)  as  the 
difference  between  the  actual  gravity  vector  and  the  reference  gravity  vector,  has  a  root 
mean  square  (RMS)  value  over  the  United  States  of  17  yg  in  the  gravity  anomaly  or  down 
component  {5).  Extremes  in  excess  of  100  yg  occur.  While  the  RMS  values  are 
representative  of  the  gravity  disturbance,  their  actual  values  are  very  position  dependent, 
varying  with  longitude  as  well  as  latitude  and  height.  Further,  their  correlation  distances 
are  rather  short,  about  20  km  for  the  east  and  north  components;  the  position  dependence 


an  ellipsoid.  In  terms  of  spherical  harmonics,  the  componets  of  jyare  given  by: 

5  J4  L  -  30  cos  L  +  3)  I 

8  “ 


Ye  =  0 

Yn  =  -  3GM  ^  ^ 

(7  cos^  L  -  3)1 


(3) 


J2  +  5  J4  R 

6  r 


where 

GM  =  gravitational  constant  times  mass  of  Earth 
L  =  geocentric  latitude  of  vehicle  (approximated  by  geodetic  latitude) 

R  =  equatorial  radius  of  earth 

32»34  =  second  and  fourth  degree  zonal  harmonic  coefficients  of  geopotential 

The  earth  rotation  term,  fl.  ft.  r,  is  computed  along  with  the  reference  gravita- 
tion  vector  so  that  we  may  consider  the  gravity  vector  as 

g=2*+A£ 

* 

where  =  reference  gravity  vector. 

The  components  of  the  reference  gravity  vector  are  given  by 

*  22 

Y  u=^u*'“  ie""®"  ^ 

Y*  =0  (5) 

'  e 

*  * 

Y  =  Y  -  tu  .  sin  L  cos  L 
'  n  'n  le 

u.  =  earth  rotation  rate 
le 

The  gravity  disturbance  vector,  6g,  which  can  be  defined  via  equation  (4)  as  the 
difference  between  the  actual  gravity  vector  and  the  reference  gravity  vector,  has  a  root 
mean  square  (RMS)  value  over  the  United  States  of  17  pg  in  the  gravity  anomaly  or  down 
component  f5).  Exfremes  in  excess  of  100  pg  occur.  While  the  RMS  values  are 
representative  of  the  gravity  disturbance,  their  actual  values  are  very  position  dependent, 
varying  with  longitude  as  well  as  latitude  and  height.  Further,  their  correlation  distances 
are  rather  short,  about  20  km  for  the  east  and  north  components;  the  position  dependence 


is  seen  to  be  rather  strong.  Thus,  it  is  advisable  to  incorporate  into  the  gravity 
computation  feedback  loop  a  computation  of  the  gravity  disturbance  as  a  function  of 
position. 


A  rigorous  computation  of  the  gravity  disturbance  components  involves  a  world 
wide  integration  of  gravity  anomaly  data  via  the  equations  (6): 


The  computer  storage  and  time  needed  to  calculate  from  these  equations  make 
their  real-time  use  impractical.  The  alternative  to  this  computation  is  the  development 
of  modeling  techniques  in  which  the  gravity  disturbance  vector  can  be  accurately 
represented  by  some  mathematical  function  which  is  adaptable  to  real-time  computation 
and  is  conservative  of  avionics  computer  resources. 

Although  several  mathematical  models  may  be  applicable  to  the  gravity 
disturbance  modeling  problem,  two  techniques  have  shown  particular  promise  in  investiga¬ 
tions  conducted  at  the  Defense  Mapping  Agency  Aerospace  Center.  The  modeling 

technique  which  has  withstood  the  test  of  time  and  application  is  the  use  of  point  masses 
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(7).  A  relatively  recent  technique,  still  under  investigation  at  DMA  Aerospace  Center,  is 
the  Finite  Element  Method  in  which  Chebyshev  polynominals  are  used  to  model  the 
gravity  disturbances  in  small  volumes  (finite  elements)  (s).  See  references  (9]  and  (lO)  for 
results  of  DMA  Aerospace  Center  investigation  into  goodness  of  fit  of  the  Finite  Element 
modeling  and  the  result  of  gravity  disturbance  modeling  on  a  selected  test  flight  track. 

In  the  Point  Mass  Model,  the  gravity  disturbances  are  given  by  (7) 


6 


Su 


=  GM 


N 

Z 

j=l 


1  30 

- fr 

P 


N 

=  GM  Z 
j=l 


(7) 
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GM 


where 


p  =  distance  from  jth  point  mass  to  exterior  point 
(^j/M)  =  ratio  of  jth  point  mass  to  Earth's  mass 
N  =  number  of  point  masses  in  set. 


In  the  Finite  Element  Model,  the  gravity  disturbances  are  given  by  (S). 
N  n  n-1 


6 


Su  =  ^  ^  ’'i'"!*  V-Z* 

n=0  1=0  )=o  ijk  ’ 


■ge=^ 


N  n  n-I 

Z  S  C^  T.(x,)  T.(x-)  T.  (xj 

e...  I  1  j  /  K  J 
n=o  1=0  j=o  ijk  ' 


(8) 


;  N  n  n-1 

gn=J:  £  Z  T;(x,)T.(xJT.(xj) 

"  n=o  i=o  j=o  "ijk  •  1  J  2  k 
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where 

N  =  highest  order  of  approximating  polynomials 

k  =  n-i-j 

Xi  =(h-h^j-)/Ah 

Xj  . 

X3  =  (L-Lr^j)ML 

T^(x)  =  shifted  Chebyshev  polynomials 

C  ,  C  ,  C  =  coefficients  of  the  set 
'^ijk  ®ijk  ijk 

^ref'  ^ref'  ^ref  =  coordinates  of  lower  southwestern  point  of  cell 

A  h,  A  X,  A  L  =  dimensions  of  the  cell 

In  both  the  Point  Mass  or  Finite  Element  techniques,  the  parameters  of  the  fit, 

(m  /M)  or  C  ,  C  ,  C  ,  are  determined  by  least  square  fit  to  available  geopotential 
^  '^ijk  ®ijk  "ijk 

data. 

In  our  previous  investigations  of  the  Finite  Element  Method  (9),  we  selected  a  9°  x 
15°  X  9  km  region  (32°N  -  41  °N,  IQ8°W  -  93°W,  0  -  9  km  altitude)  in  the  southwestern 
part  of  the  continental  United  States  (Figure  2)  to  perform  test  computations.  It  was 
found  that  third  order  polynomials  would  suffice  to  model  the  east  and  north  components 
of  the  gravity  disturbance  vector  to  an  accuracy  of  +  5  yg  in  102  elements  in  this  area. 
Of  the  remaining  elements,  29  were  adequately  modeled  with  sixth  order  polynomials 
(stippled  on  Figure  2)  while  the  remaining  four  required  ninth  order  polynomials  (stippled 
and  boxed  on  Figure  2).  Along  a  simulated  flight  line  which  carries  an  aircraft  at  Mach 
0.5  over  a  mixture  of  cells  requiring  third,  sixth  and  ninth  order  polynomials,  we  found 
that  if  third,  sixth  and  ninth  order  polynomials  were  used  as  required,  the  modeling 
accuracy  achieved  was  +  2.4  y  g  in  the  east  component  and  +2.7  y  g  in  the  north  component. 
This  accuracy  will  be  accepted  in  this  paper  as  representative  of  that  which  can  be 
achieved  by  either  Point  Mass  or  Finite  Element  modeling. 


III.  ERROR  PROPAGATIONS  AND  KALMAN  FILTERING 

A  semi-continucxis  Kalman  Filter  process  is  used  in  which  the  error  state  vector  is 
propagated  by  the  equation 


Where 


X  =  F(t)  X  +  W  +  I! 


X  =  error  state  vector 


F(t)  =  fundamental  matrix,  essentail^y  Fx 
W  =  vector  of  error  forcing  functions 
U  =  vector  of  correcting  functions 

Assuming  that  the  vertical  channel  is  damped  by  a  barometric  or  radar  altimeter,  the 
state  vector  consists  of  east  and  north  position  errors,  east  and  north  velocity  errors  and 
platform  orientation  errors  about  the  east,  north  and  down  axes. 

The  error  covariance  matrix  is  progagated  by  the  equation 


P  =  F(t)P  +PF‘(t)  +  Q(t) 


where 


P  =  error  covariance  matrix  of  system 

Q(t)  =  error  covariance  of  forcing  functions  minus 
correcting  functions  if  any 

Note  that  W  is  considered  to  be  white  noise  with  zero  mean  and  covariance  Q(t). 

A  discrete  Kalman  Filter  is  introduced  to  incorporate  position  and  velocity 
updates,  data  which  will  be  available  at  discrete  intervals.  The  error  state  vector  is 
updated  by  the  equation 

X^"^^  =  -  K(Z  -  (1 1) 

and  the  error  covariance  matrix  by  the  equation 

P^^^  =  (l- KH)P^“^  (12) 


where 


Z  =  a  measurement  matrix 


H  =  a  matrix  relating  the  measurement  to  the  state  vector 
in  the  sense  Z  =  HX 
K  =  the  Kalman  gain  matrix 

The  matrix  Z  is  assumed  zero,  since  the  measurement  resets  the  position  and  velocity 
with  "zero"  error. 

The  Kalman  gain  matrix  is  the  keystone  of  the  filtering  process  in  that  it  provides  for 
relative  weighting  between  the  error  covariance  of  the  state  vector  and  the  error 
covariance  of  the  measurement.  It  is  given  by: 

K  =  -  R] 

where 

R  =  error  covariance  matrix  of  measurement. 

Note:  In  equations  (II)  through  (13),  the  superscript  implies  "immediately  before 

update",  and  the  superscript  implies  "immediately  after  update." 

A  gyro  drift  rate  of  0.00000 l°/hour  and  an  accelerometer  bias  of  1  pgwere 

included.  In  each  error  propagation,  the  vector  of  error  forcing  functions,  W,  includes  not 

only  these  effects  but  also  the  unmodeled  gravity  disturbance  vector  6  g,  which  we  have 

computed  for  all  points  on  the  simulated  flight  line  by  application  of  equations  (6)  to 

actual  gravity  data.  The  error  covariance  Q  of  this  forcing  function  is  a  diagonal  matrix 

10  2 

with  the  gravity  disturbance  variances  set  to  64  x  10"  (m/sec/sec)  for  the  unmodeled 
gravity  disturbances  and  equal  to  the  modeling  accuracy  stated  above  for  the  modeled 
gravity  disturbances  (ll). 

Two  propagations  were  accomplished  over  the  flight  line  depicted  on  Figure  2  to 
establish  a  baseline  reflecting  the  effect  of  the  gravity  disturbances  with  no  updating  by 
GPS  data.  One  propagation  was  accomplished  with  the  vector  of  correcting  functions,  U, 
equal  to  zero,  representing  the  error  state  vector  which  would  result  if  the  gravity 
disturbances  were  left  unmodeled.  The  second  propagation  was  accomplished  with  U  set 
equal  to  the  value  determined  by  the  Finite  Element  method  as  a  function  of  position. 
Results  of  these  propagations  are  given  In  Table  I. 


8 


Component 

TABLE  1 

Maximum  Error  State  Vector  (no  GPS  Update) 

Unmodeled 

Modeled 

E 

1059.24m 

31.78m 

N 

1450.59m 

62.10m 

• 

E 

1.218m/sec 

0.0  38  m/sec 

N 

1.598m/sec 

0.069m/sec 

Three  additional  sets  of  propagations  were  accomplished  in  order  to  assess  the 
impact  of  the  GPS  data  on  the  inertial  navigation  system,  with  and  without  the  gravity 
disturbance  modeling.  In  these  propagations  the  measurement  covariance  matrix  R  is  a 
diagonal  matrix  with  the  pertinent  elements  set  equal  to  the  accuracies  cited  earlier  in 
this  paper  for  the  GPS.  A  measurement  is  assumed  available  at  each  time  increment 
along  the  flight  line.  The  inertial  navigator  may  use  the  position  data,  the  velocity  data, 
or  both  for  updating.  Table  D  gives  the  results  for  the  case  where  both  are  used. 

TABLE  II 


Maximum  Error  State  Vector  (Position  and  Velocity  Data) 


Component 

Unmodeled 

Modeled 

E 

73.97m 

8.08m 

N 

42.43m 

8.44m 

• 

E 

0.233m/sec 

0.017m/sec 

N 

0.162m/sec 

0.019m/sec 

The  application  of  the  GPS  data  has  improved  the  unmodeled  case  by  better  than  an  order 
of  magnitude  in  position  and  nearly  as  well  in  velocity.  However,  the  improvement  is  not 
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as  great  as  that  afforded  by  gravity  disturbance  modeling  without  the  GPS  data  (see  Table 
I).  When  the  inertial  navigator  is  furnished  with  both  GPS  data  and  gravity  disturbance 
modeling,  the  results  are  quite  dramatic. 


Tables  III  ahd  IV  give  the  results  for  the  cases  where  position  data  only  and 
velocity  data  only,  respectively,  are  used. 

TABLE  III 


Component 

Maximum  Error  State  Vector  (Position  Data  Only) 

Unmodeled 

Modeled 

E 

82.32m 

8.46m 

N 

47.90m 

8.89m 

• 

E 

0.255m/sec 

0.017m/sec 

N 

0. 177m/sec 

0.019m/sec 

Component 

TABLE  IV 

Maximum  Error  State  Vector  (Velocity  Data  Only) 

Unmodeled 

Modeled 

E 

248.28  m 

21 .88m 

N 

220.92m 

20.58m 

« 

E 

0.284m/sec 

0.0  28  m/sec 

N 

0.236m/sec 

0.026m/sec 

The  use  of  position  information  alone  to  update  the  inertial  navigator  yields  error  state 
vectors  nearly  the  same  as  those  using  both  position  and  velocity  data.  The  velocity 
information  does  not  contribute  as  significantly  to  the  updating  process;  however,  velocity 
data  provides  a  substantial  improvement  in  the  case  where  gravity  disturbances  remain 
unmodeled. 

Figures  3  through  6  graphically  portray  the  variations  of  the  error  as  a  quasi- 
time-dependent  accelerometer  bias.  On  these  graphics,  which  are  computer  drawn,  the 
fainter  lines  are  the  position  and  velocity  errors  due  to  the  unmodeled  gravity  disturbance 


vector,  while  the  heavier  lines  are  the  position  and  velocity  errors  after  compensation  for 
the  gravity  disturbances  has  been  applied. 

Because  the  updating  occurs  in  a  semi-continuous  mode,  i.e.,  at  each  time 
increment  along  the  flight  line,  we  do  not  see  glitches  or  a  saw-tooth  effect  in  the 
graphics.  Instead,  the  graphics  approximate  sine  waves  with  a  frequency  equal  to  the 
Schuler  frequency  (period  =  minutes).  Higher  frequency  variations,  especially  apparent 
in  the  velocity  graphic,  reflect  the  effect  of  the  short  wavelength  components  of  the 
gravity  field. 

IV.  CONCLUSIONS 

Although  the  GPS  provides  the  most  accurate  and  effective  means  of  providing 
updated  position  and  velocity  data  to  an  inertial  navigation  system,  the  error  state  vector 
may  grow  to  greater  limits  than  that  tolerable  in  certain  applications.  Using  optimum 
gravity  disturbance  modeling  techniques,  without  further  updating  by  the  GPS,  the 
maximum  error  state  vector  can  be  made  comparable  with  that  obtainable  by  the  GPS.  It 
is  seen  then  as  an  alternative  should  GPS  be  denied  to  the  inertial  navigator.  The  fullest 
capabilities  of  both  are  realized  when  the  inertial  navigator  has  both  gravity  disturbance 
modeling  and  position  and  velocity  updating  from  the  GPS.  If,  in  the  interest  of 
conserving  computer  resources,  it  is  desirable  to  use  only  position  data  in  the  updates, 
there  is  only  a  minor  degradation  in  performance  of  the  system.  Use  of  only  velocity  data 
in  the  updates  offers  an  improvement  over  no  modeling  at  all. 
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